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Preface 


This  final  report  concerning  modeling  semiconductor  and  quartz  crystal 
growth  is  focused  on  the  transport  equation  of  the  dopant  and  of  the  nutrient 
in  neighborhood  of  the  melt/solid  and  liquid/solid  interface. 

We  present  the  weakness  of  the  "stagnant-film"  or  "diffusion-layer"  model 
and  using  recent  developments  in  the  field  of  periodic  material  modeling  we 
identify  new  transport  equation  in  the  region  adjacent  to  the  interface.  This 
equation  takes  into  account  the  interaction  between  complicated  global  flow 
patterns  and  the  convective-diffusive  dispersion  mechanism,  which  exist  in 
the  so-called  “mushy  zone”  or  “precrystallization  region”.  We  assume  that  in 
this  zone  we  have  a  periodic  structure,  which  is  a  weak  form  of  the  periodic 
structure  of  the  crystal.  If  the  periodic  structure  in  the  precrystallization 
region  is  no  significant  i.e.  the  size  of  the  periodically  distributed  solid 
inclusions,  which  are  small  with  respect  to  the  distance  between  two 
neighbouring  inclusions,  is  smaller  than  a  critical  size;  then  global  flow  is 
not  influenced  and  the  transport  equation  of  the  solute  is  the  classical 
convective-diffusive  equation  except  the  diffusive  term  which  is  more 
refined  because  a  diffusion  tensor  appears.  But  if  the  periodic  structure  in  the 
precrystallization  region  is  significant  i.e.  the  size  of  the  periodically 
distributed  inclusions  is  greater  than  a  critical  size;  then  the  flow  in  mushy 
zone  is  a  Darcy  flow  and  in  the  new  transport  equation  of  the  solute  we  have 
new  more  refined  convective,  diffusive  terms.  It  must  be  noted  that  these 
new  equations  do  not  take  into  account  the  interaction  of  the  dopant  field 
with  the  crystal.  A  kind  of  interaction  appears  just  in  boundary  conditions.  In 
this  context  we  can  think  to  improve  our  equations  by  the  addition  of  a  new 
term,  which  expresses  this  interaction. 

Concerning  the  accuracy  of  our  equations  we  want  to  point  out  that  we  have 
established  a  convergence  result,  which  provejthat  our  equations  obtained  by 
homogenization  are  an  approximation  of  the  equations  governing  the  real 
phenomena.  But,  in  this  context  it  is  obviously  that  in  the  future  it  will  be 
necessary  to  have  also  numerical  evidences. 


Timisoara  Stefan  Balint 

December  1998 


1.  On  the  determination  of  the  radial  segregation  of  dopant  in 
the  case  of  semiconductor  crystal  growth  in  Bridgman  - 
Stockbarger  configuration 


1.1  Introduction 


The  compositional  uniformity  of  crystals  grown  from  the  melt  depends 
strongly  on  the  pattern  and  intensity  of  flow  in  the  melt  and  on  the  shape  of 
the  melt/solid  interface.  In  the  processes  such  as  the  Czochralski  and  floating 
zone  methods  and  the  vertical  Bridgman- Stockbarger  configuration  studied 
here,  the  unequal  partitioning  of  dopant  between  melt  and  crystal  during 
solidification  causes  a  concentration  gradient  perpendicular  to  the  melt/solid 
interface.  When  this  interface  is  planar  and  the  melt  is  quiescent,  except  for 
motion  caused  by  crystal  growth,  the  concentration  field  decays 
exponentially  with  distance  into  the  melt.  Curvature  of  the  solidification 
interface  and  convection  in  the  melt,  both  change  the  concentration 
distribution  along  the  interface  and  alter  the  dopant  level  in  the  crystal. 

The  influence  of  convection  on  segregation  in  the  growth  direction  (axial 
segregation)  was  pointed  out  by  the  Burton’s  et  al.  analysis  [1-3]  of  the 
effect  of  crystal  rotation  on  mass  transfer.  When  the  crystal  diameter  is  large 
and  the  velocity  field  is  similar  along  the  radius  of  the  melt/solid  interface, 
convection  only  alters  the  concentration  field  perpendicular  to  the  interface. 
Burton  et  al.  analyzed  this  case  and  developed  an  expression  for  an  effective 
segregation  coefficient  and  introduced  the  notion  of  an  axial  “boundary  — 
layer”  thickness  for  diffusion-controlled  mass  transfer.  Others 
[4-8]  applied  the  idea  of  a  uniform  diffusion-layer  or  stagnant-film  that  is 
masking  a  growing  crystal  from  a  well  mixed  bulk  and  have  determined  so- 
called  boundary-layer  (more  appropriately  diffusion-layer)  thickness  without 
any  picture  for  the  fluid  motion  in  the  melt.  Other  than  Burton,  Prim  and 
Slichter’s  original  analysis  of  rotating  flow  few  papers  dealing  with  melt 
crystal  growth  [9,10]  have  addressed  the  exact  coupling  between  the  fluid 
flow  and  mass  transfer  in  a  consistent  way. 

In  many  growth  configurations,  the  concept  of  a  uniform  diffusion-layer 
yields  an  over-simplified  view  of  the  role  of  convection  in  dopant 
segregation  [11].  When  the  flow  is  laminar  and  cellular,  as  in  the  case  for 
many  small-scale  growth  systems  and  for  reduced  gravity  experiments. 
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convective  mass  transfer  is  uneven  along  the  surface  of  the  crystal  and 
significant  lateral  segregation  results.  The  coupling  between  moderate 
convection  and  lateral  dopant  segregation  thoroughly  analyzed  for  the 
rotational ly-driven  flows  in  small-scale  floating  zones  [10,12,13]  and  has 
been  demonstrated  for  buoyancy-drive  convection  under  microgravity 
conditions  [14].  All  this  studies  were  of  model  crystal  growth  systems  with 
planar  solidification  interfaces. 

Coriell  and  Sekerka  [15]  (also  see  ref.  [16])  demonstrated  that  significant 
radial  segregation  occurs  in  systems  without  convection  tangent  to  the 
crystal  surface  when  the  radius  of  curvature  of  this  interface  is  the  same  or 
less  than  the  length  scale  of  the  concentration  gradient  adjacent  to  the 
interface.  Curvature  induced  segregation  was  shown  to  be  an  important 
contribution  to  dopant  inhomogeneity  in  capillary  growth  systems  [17-19] 
and  may  be  important  in  other  small-scale  crystal  growth  experiments  in 
which  convection  in  the  melt  has  been  suppressed. 

We  consider  the  prototype  vertical  Bridgman  -  Stockbarger  growth  system 
presented  in  Fig.l. 
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Fig.l.  Geometry  of  the  prototype  vertical  B-S  system 
The  influence  of  the  natural  convection  governed  by  the  equations 
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(1)  — ^  + v-Vv  =  -Vp  +  Pr-V^  v-Ra-Pr-0-e3 
ot 

(2)  — +  v-ve  =  v^0 

a 

(3)  V-v  =  0 

on  the  dopant  dispersion  governed  by  the  equation 

Wdc  Sc  -  _ 
a  Pr 

was  analyzed  by  C.J.  Chang  and  R.A.  Brown  in  [20].  The  calculations 
described  here  are  of  the  axisymmetric  steady  state  velocity,  temperature  and 
concentration  fields  and  melt/solid  interface  shapes.  The  mathematical  free¬ 
boundary  problem  that  describes  these  variables  is  solved  by  finite- 
element/Newton  technique  that  computes  simultaneously  the  shape  of  the 
solidification  interface,  the  velocity  and  pressure  fields  in  the  melt  and  the 
temperature  distribution  in  both  melt  and  crystal.  Flow  field  calculated  with 
the  finite-element/Newton  algorithm  are  used  in  a  separate  calculation  of  the 
distribution  of  a  dilute  dopant  within  the  melt  and  crystal.  The  bulk  of  the 
calculations  reported  here  are  for  a  melt  and  crystal  with  thermophysical 
properties  similar  to  those  of  the  gallium-doped  germanium  system. 

Three  distinct  types  of  flow  patterns  were  observed.  At  low  Rayleigh 
numbers  (Ra  <  10^)  the  streamlines  were  rectilinear  and  only  slightly 
distorted  by  buoyancy  forces.  For  intermediate  values  of  Ra  (10^  <  Ra  <  10^) 
a  cellular  flow  developed  which  was  driven  by  the  radial  temperature 
gradients  established  by  the  mismatch  in  thermal  boundary  condition^ 
between  the  adiabatic  and  hot  zones.  The  flow  moved  upward  along  the 
sidewall  and  downward  at  the  centerline  of  the  ampoule.  The  center  of  the 
cell  was  located  slightly  above  the  gradient  zone  and  migrated  downward 
and  toward  the  sidewall  with  increasing  Ra.  Increasing  the  Rayleigh  number 
to  2.6  X  10^  led  to  the  development  of  a  weak  secondary  cell  adjacent  to  the 
melt/solid  interface.  The  outset  of  the  multi-cellular  flows  marks  the  start  of 
the  third  type  of  flow  pattern.  The  motion  in  the  secondary  cell  next  to  the 
interface  was  in  the  opposite  direction  to  the  main  cell  and  led  to 
qualitatively  different  radial  segregation  than  the  flows  with  a  single  cell 
which  exist  for  Ra  less  than  1x10^. 
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The  temperature  fields  show  the  effect  of  the  Prandtl  number  in  the  gallium 
germanium  system.  For  Rayleigh  numbers  between  zero  and  1  x  lO'*  the 
thermal  fields  in  both  the  melt  and  crystal  were  essentially  the  same  as  the 
field  calculated  without  convection  (Ra=0)  and  were  similar  to  results  of  Fu 
and  Wilcox  [21].  Increasing  Ra  above  1  x  10'*  caused  the  isotherms  along 
the  axis  of  the  melt  to  compress  toward  the  melt/solid  interface  by  the 
downward  fluid  motion.  By  Ra  =  5  x  10^  the  shape  of  several  isotherms 
farthest  from  the  interface  inverted  from  convex  to  concave  at  the  center  of 
the  melt.  The  isotherms  in  the  crystal  were  unchanged  by  changing  Ra.  Also, 
the  fact  that  a  large  position  of  the  crystal  was  at  the  uniform  temperature 
0=0. 1  demonstrated  that  the  length  of  ampoule  in  the  cold  portion  of  the 
furnace  was  sufficient  to  guarantee  that  the  position  of  the  end  of  the 
ampoule  was  unimportant. 

The  shape  of  the  melt/solid  interface  was  unchanged  by  convection  for  Ra 
between  0  and  10^  For  higher  Rayleigh  numbers,  the  hot  melt  moving  down 
the  axis  of  the  ampoule  drove  the  melt/solid  interface  deeper  into  the 
adiabatic  region.  The  changes  in  interface  shape  caused  by  convection  were 
not  large;  even  for  Ra  =  5  x  10^  the  deflection  of  the  interface  was  only  6% 
of  its  mean  location. 

Dopant  fields  were  calculated  for  the  velocity  field  and  melt/solid  interface 
shapes  discussed  above  for  the  segregation  coefficient  k=0.1  and  Schmidt 
number  Sc=10  similar  to  the  gallium-doped  germanium  system.  The  almost 
parallel  iso-concentration  lines  for  Rayleigh  numbers  up  to  10  correspond  to 
the  one-dimensional  solidification  model.  The  concentration  field  was 
deformed  at  higher  values  of  Ra;  at  Ra  =  1  x  10'*  the  concentration  field  had 
the  beginnings  of  the  uniform  core  of  melt  with  steeps  concentration 
gradients  along  each  boundary,  consistent  with  the  boundary-layer  model  for 
a  well-mixed  melt.  A  large  amount  of  radial  segregation  was  present  even  at 
this  level  of  convection:  flow  downward  along  the  axis  swept  dopant  from- 
crystal  and  induced  radial  segregation  across  the  surface  of  the  crystal.  The 
change  in  concentration  across  the  interface  evolved  with  increasing  Ra 
from  approximately  1%  segregation  caused  by  interface  curvature  at  Ra=0  to 
almost  70%  of  the  mean  value  1/k  for  the  worst  case  of  Ra  =  1  x  10^  The 
radial  variation  of  dopant  decreased  with  the  more  intense  flow  motion  that 
corresponded  to  Ra=10'^.  The  maximum  in  radial  segregation  with  increasing 
convection  found  here  was  consistent  with  calculations  of  Nikitan  for  the 
growth  of  gallium-doped  germanium  in  a  horizontal  boat. 

The  percent  radial  segregation,  defined  as  Ac=[c(0,h(0))-c(A,h(A))]xl00/k 
reached  a  maximum  for  the  flow  comesponding  to  nearly  Ra=10^  and 
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decreased  for  larger  values  of  Rayleigh  number.  The  value  of  the 
segregation  coefficient  set  the  mean  concentration  in  the  melt  at  the  interface 
and  affected  the  level  of  radial  segregation  in  the  crystal.  The  exponentially 
decreasing  concentration  profile  present  at  low  values  of  Ra  extended  further 
into  the  melt  for  smaller  segregation  coefficients  and  made  the  concentration 
field  more  sensitive  to  convection.  Consequently  the  level  of  segregation  Ac 
was  higher  for  lower  k. 

Changing  the  diffusivity  of  the  solute  also  had  a  marked  effect  on  the  level 
of  radial  segregation.  At  high  Schmidt  numbers  and  Rayleigh  numbers  of 
1x10  and  above,  the  isoconcentration  curves  developed  fingers  oriented 
parallel  to  the  melt  solid  interface  that  corresponded  to  extremely  rapid 
variation  in  composition  within  a  distance  of  the  same  order  of  magnitude  at 
the  thickness  led  to  large  (over  1 00%)  radial  segregation. 

The  cornplicated  segregation  patterns  discussed  created  doubts  as  to  the 
applicability  of  the  simple  “diffusion-layer”  or  stagnant-film  model  for 
correlating  axial  segregation  behavior.  Details  results  for  the  velocities  and 
dopant  profile  make  possible  to  check  this  correlation.  To  do  this,  radially 

A 

averaged  concentration  profiles,  defined  as  c(z)  =  Jc(r,  z)  •  r  dr  was 

0 

calculated  for  different  values  of  Ra,  k  and  Sc.  For  law  Rayleigh  numbers 
(0<  Ra  <  1  X  10^),  these  radially  averaged  profiles  were  essentially  the  same 
as  the  concentration  profiles  predicted  by  one-dimensional  models  which 
account  only  for  the  convection  caused  by  crystal  growth.  A  region  of  nearly 
uniform  average  concentration  c(z)  =  Cg  developed  with  increasing  Ra.  The 

profile  for  Ra=l  x  10“^  was  divided  into  a  region  of  thickness  5f  with  steep 
dopant  gradient  adjacent  to  the  crystal,  a  zone  of  uniform  concentration  and 
a  gradient  zone  caused  by  the  fictitious  inlet  condition  at  the  top.  These_ 
profiles  were  reminiscent  of  the  form  assumed  by  a  stagnant  film  mode^and  ^ 
suggest  the  equivalence  of  5f  with  the  width  of  usual  diffusion  layer.  The 
thickness  5f  were  compared  with  results  of  this  approach  through  calculation 
of  effective  segregation  coefficient.  The  effective  segregation  coefficient 
appropriate  for  a  true  unsteady  Bridgman  system  with  A=0.25  was 
appioximated  as  keff=Cs/cB=l/cB  where  Cs  was  the  dimensionless  average 
concentration  across  the  crystal. 

If  a  diffusion  layer  of  thickness  5  =  5/L  existed  that  separated  the  crystal 
surface  from  bulk  melt  at  concentration  Cq,  the  effective  segregation 
coefficient  keff  would  be  derived  following  Flemings  formula  as 
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1  —  s  _ 


k  +  (l  -  k;)exp(-  Vl  •  6  /  d) 
k 


or  in  dimension  less  form  used  here 


as  kgff  !<;  + (i-k)exp(-a-5-Pe-Sc/Pr) 

Rearranging  this  yielded  an  explicit  relation  for  the  diffusion-layer  thickness 
/  ,  ,  1.  ^ 


5  =  - 


Pr 


In 


1-k 


eff 


.  Values  of  5  computed  from  the  above 

Pe-Sc-G  J 

formula  with  kgff  calculated  from  the  finite-element  results  are  compared  to 
values  of  5f  measured  from  the  profiles  of  c(z).  In  each  case,  5  and  5c  were 
of  the  same  order  of  magnitude  but  different  by  as  much  as  a  factor  of  two. 
This  result  is  not  surprising.  Within  the  distance  5r  from  the  crystal  in  the 
finite  element  calculation,  the  velocity  field  was  a  combination  of  the  growth 
velocity  and  a  contribution  from  natural  convection.  These  two  components 
were  of  the  same  order  of  magnitude  and  it  was  the  latter  component  which 
caused  the  discrepancy  between  5  and  5f  which  led  to  radial  segregation. 

One  of  the  main  conclusions  of  Chang  and  Brown  in  the  paper  [20] 
concerning  the  steady  state  is  the  following: 

“For  moderate  levels  of  convection  the  “stagnant  film”  or  “diffusion-layer” 
model  is  a  gross  oversimplification  of  the  interaction  between  complicated 
flow  patterns  and  the  dopant  field.  Diffusion-layer  thicknesses  5  determined 
by  the  formula 


Pr 


•In 


k  1-k 


eft 


1-k 


(5)  5  = 

Pe  •  Sc  •  G  k*^eft  i  J 
and  experimental  data  are,  at  the  best,  empirical  fits  to  effective  segregation 
coefficient,  for  the  crystal/melt  system.  The  comparison  demonstrates  that, 
although  the  radial  averaged  diffusion-layer  thickness  determined  fiom  the 
finite-element  simulations  and  formula  (5)  are  of  the  same  order  of 
magnitude,  the  actual  concentration  gradient  adjacent  to  the  crystal  is  far 
from  radial  uniform.  As  much  60%  radial  segregation  can  exist.  Only 
detailed  calculations  of  the  exact  interaction  of  fluid  flow  and  dopant  profile 
adjacent  to  the  interface  can  be  used  to  estimate  the  level  of  radial 
segregation  in  crystal.  Our  study  of  a  prototype  Bridgman  system  gives 
qualitative  understanding  of  fluid  flow  and  dopant  segregation  in  actual 
growth  system  and  will  serve  as  the  starting  point  for  more  refined 
calculations,  aimed  at  comparison  with  experiment  in  well-characterized 
small-scale  system”. 
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Starting  from  this,  we  propose  a  new  model  in  which  the 
"precrystallization  zone"  or  "mushy  zone"  substitutes  the  "stagnant 
film"  or  "diffusion-layer"  like  M.B.  Tahar  in  [22].  Tahar  approach  this 
solid-liquid  mushy  zone  using  a  model  of  a  porous  media  with  evolving 
heterogeneity,  and  adapt  a  spatial  averaging  method  to  derive 
macroscopic  equations  describing  geometry  of  the  dendritic  structure  as 
well  as  dispersive  effects.  For  as  the  “mushy  zone”  is  a  thin  bed  of 
periodically  distributed  solid  inclusions,  which  are  small  with  respect  to 
the  distance  behveen  two  neighboring  inclusions  (we  shall  say  that  the 
concentration  is  small). 

The  diffusion-layer  thicknesses  5  are  substituted  by  the  “mushy  zone” 
thicknesses  which  can  be  determined  by  X-ray  diffraction 
measurements.  In  this  way  the  diffusion-layer  or  stagnant  film  masking 
the  growing  crystal  is  substituted  by  a  thin  periodic  porous  medium 

bounded  below  by  an  impermeable  rigid  wall  which  is  the  melt/solid 
interface. 

If  the  size  of  solid  inclusions  is  smaller  like  the  critical  size  (see  [23])  the 
flow  in  “mushy-zone”  is  the  global  flow  given  by  (l)-(3). 

If  the  size  of  solid  inclusions  is  critical  (see  [23])  the  flow  in  the  mushy 
zone  is  a  Brinkman  flow,  determined  by  the  global  pressure  and 
respecting  at  the  “mushy-zone”  and  “pure-fluid”  interface  the  Beavers- 

Joseph  or  Jones  modified  Beavers-Joseph  or  Rudriah  conditions  tsee 
[24],[25],[26]).  ^ 

If  the  size  of  solid  inclusions  is  larger  like  the  critical  size  (see  [23])  then 
the  flow  in  the  “mushy-zone”  is  a  Darcy  flow,  determined  by  the  global 
pressure  and  respecting  at  the  “mushy-zone”  and  “pure-fluid”  interface 
Rudriah  conditions  (see  [26]). 

The  dispersion  mechanism  of  the  dopant  in  mushy-zone  is  due  to  the 
diffusion  and  convection,  which  exists  simultaneously,  and  is  assumed 
that  balance  each  other  or  convection  dominates  diffusion.  Obviously 
the  mechanism  takes  into  account  also  the  existence  of  some  solid 
inclusions.  This  dispersion  mechanism  model  for  the  dopant  adjacent  to 
the  interface  is  more  refined  as  the  dispersion  mechanism  in  diffusion- 
layer  model  as  well  as  the  dispersion  mechanism  described  by  equation 

(4)  and  is  governed  by  a  new  convective-diffusive  equation  which 
generalize  equation  (4). 

In  the  new  convective-diffusive  equation  we  use  the  velocity  field  the 
melt/solid  interface  defined  by  (l)-(3)  and  a  certain  periodic 
microstructure  similar  with  the  microstructure  of  the  crystal. 
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We  find  also  the  stationary  solutions  without  radial  segregation  of  the 
new  convective-diffusive  equation  and  perturbations  with  respect  to 
these  solutions  are  asymptotically  stables.  We  establish  also  the  optimal 
control  for  the  stationary  solutions  without  radial  segregation. 


1.2  The  flow  in  the  mushy  zone 

For  a  periodic  material  QczR^  for  which  the  period  is  a  cubic  cell  homotetic 
with  a  small  ratio  S  of  the  unit  cell  Y  =  [-1/2,  1/2]'cr\  in  which  the  fluid 
domain  Yp  and  the  solid  on  Ys  have  a  smooth  boundary  F  (like  in  Fig.2.), 
using  homogenization  method  it  was  deduced  a  macroscopic  model,  from 
microscopic  processes,  assuming  that  the  flow  is  slow  and  the  size  of  the 
solid  inclusions  is  of  the  same  order  like  the  size  of  the  cell  (low  porosities) 
(see  [27],  [28],  [29]). 


Ys 


Yf 


Fig.2.  The  unit  cell  Y 

In  the  macroscopic  model  the  flow  <  v  >  (x)  is  given  by  Darcy  s  law 
(6)  <v>(x)  =  -K*  VxP®  ,  xeQ 

where  K  is  the  permeability  tensor  defined  by  the  coefficients  Ky  given  by 

and  v^  are  the  components  of  the  micro  velocity  fields  defined  in  the  unit 
cell  Y  by  the  following  boundary  value  problem. 


9 


j  =  1,2,3 


jiAyvJ  =  VyqJ  -  e:  in  Yj 


Vyvj=0  inYp 

V  j  =0  on  r 

v-^  and  qJ  are  Y  periodic 
Cj  is  the  unit  vector  in  yj 


direction 


In  formula  (6)  represents  the  pressure  gradient  in  Q,  Numerical 
calculus  of  the  permeability  tensor  coefficients  Ky,  solving  (8)  by  finite- 
elements  method,  was  performed  by  Chiang.  C.  Mei,  Jean-Luis  Auriault  and 
Chiou-On  Ng  in  the  paper  [30]  and  by  Cheo  K.  Lee,  Chin-Cheng  Sun  and 
Chiang  C.  Mei  in  [31].  They  consider  the  case  when  the  solid  inclusion  is  a 
Wigner-Seitz  grain  in  the  unit  cell  Y.  In  this  case  due  to  the  symmetry, 
Kij  =  0  for  i  ^  j  and  Kn  =  K22  =  Kss-  Therefore  it  is  sufficient  to  consider 
only  Kn  induced  by  the  unit  pressure  gradient  in  e,  direction.  Computation 

Iy  I  |y  I 

have  been  performed  for  porosities  ttv  =  in  the  ranee 

|Y|  |Y| 

0.2518  <  Tty  <  0.8333  and  the  results  were  compared  to  the  values  of  the 
permeability  given  by  the  well  known  empirical  Kozeny-Carman  formula 


(9)  k  =  --- 
5  (1 


71, 


V. 


A,  -i 


where  Vs  is  the  volume.  As  is  the  area  ofthe  grain,  and 
cube. 


£  is  one  side  of  the- 


Within  the  range  0.37  <  Tiy  <  0.68  for  which  the  empirical  formula  is  based 
on  experiments,  the  results  are  consistent  and  in  trend  with,  but  fall  slightly 
below  the  empirical  formula.  Outside  this  range  of  porosities,  the  deviation 
increases. 

It  should  be  noted  that,  the  Kozeny-Carman  formula  is  the  best  fit  to 
experimental  data  for  all  grain  shapes,  but  outside  of  the  range  of  porosities 

0.37  <%Y<  0.68  is  an  extrapolation  of  measured  data  and  may  not  be  totally 
accurate. 

In  the  paper  [31]  results  were  compared  to  numerical  values  obtained  by 
Zick  and  Hornsey  in  [32]  for  uniform  spheres  of  various  packings. 
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Discrepancies  between  the  packed  spheres  and  Wigner-Seitz  grains  appear 
for  close  packing  (low  porosity),  since  particle  interaction  is  affected  by  the 
geometry  most  significantly.  At  higher  porosity,  the  results  agree  remarkably 
well. 

Other  results  concerning  the  permeability  tensor  coefficients  in  the  case  of 
dendritic  growth  are  given  in  [22]  and  in  general  in  the  case  of 
solidifications  in  [33-45]. 

In  the  paper  [23]  Th.  Levy  studies  slow  flow  of  an  incompressible  viscous 
fluid  in  an  array  of  a  great  number  of  small  fixed  particles  assuming  that  the 

particle  size  S  and  the  distance  rj  between  two  neighboring  solids  (p  is  also 
the  size  of  the  period)  are  such  that  s  «  ri  «  1.  This  assumption 
corresponds  to  low  concentration  of  solid  inclusions  (high  porosity).  She 
proves  that  there  exists  a  critical  size  for  particles,  for  which  Brinkman’s 
flow  occurs.  For  smaller  sizes  solid  do  not  influence  the  flow  and  for  larger 
particle  sizes  the  fluid  flow  is  governed  by  Darcy's  law. 

In  this  case  it  is  worthwhile  obtaining  approximate  homogenization  formulas 
simpler  than  the  usual  ones.  This  way  be  done  by  disregarding  the 
interaction  between  inclusions:  the  local  solutions  appearing  in  the 
homogenization  formulas  (which  are  periodic,  the  periodicity  taking  into 
account  the  interaction)  become  solutions  in  the  whole  space  (without 
interaction).  In  this  connection,  Einstein  in  1906  gave  a  celebrated  formula 
for  the  homogenized  viscosity  of  a  dilute  suspension  of  rigid  solid  spherical 
particles  in  a  viscous  fluid.  As  the  general  case  of  non-small  concentration 
involves  new  phenomena  such  as  variable  microstructure,  anisotropy  and 
modification  of  the  inertial  terms,  it  is  worthwhile  obtaining  the  Einstein’s 
approximation  as  the  asymptotic  behavior  for  small  concentration  of  the 
homogenization  formulas  (see  [46],  [47]). 

Returning  to  the  mushy  zone  where  we  have  supposed  that  the  solid 
inclusions  are  small  with  respect  to  the  distance  between  two  neighboring- 
inclusions  according  to  Levy’s  paper  if  the  size  of  the  solid  particles  is 
smaller  than  the  critical  size  then  the  particles  do  not  influence  the  global 
flow,  and  therefore  in  mushy  zone  the  flow  is  the  global  flow  defined  by  the 
equations 

(10)  V-  V  V  =  -Vp°  -h  Pr  V-  Ra  •  Pr-  0  •  e'j 

(11)  v-ve  =  v-e 

(12)  Vv  =  0 
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where  V  =  er-  — +  e3-^—  is  the  gradient  operator  in  cylindrical 

coordinates.  The  Prandtl  number  Pr  and  the  Rayleigh  number  Ra  appearing 
in  equation  ( 1 0)  and  ( 1 1 )  are  defined  by 

Pr  =  -y-,Ra  =  p.g(T,-Tj.  — -H, 

at 

This  case  coiresponds  to  the  case  considered  in  [20]  and  other  papers. 

If  the  size  of  the  solid  inclusions  is  critical  then  the  flow  in  the  mushy  zone 
is  a  Brinkman  flow  governed  by  a  Brinkman’s  law 

(13)  pA<  v>=  Vp° +|Li-(p-H-<  v> 

where  cp  is  the  volume  concentration  of  solid  particles  and  H  is  the 
translation  tensor. 

If  the  size  of  solid  particles  is  greater  than  the  critical  size  then  the  flow  is  a 
Darcy  flow  governed  by  a  Darcy’s  law 

(14)  p-(p-H-<  v>=-Vp° 

The  permeability  tensor  K,  that  appears  in  critical  and  supercritical  case,  is 
expressed  in  terms  of  the  particles  concentrations  (p  and  the  translation 

tensor  H:  K  =  —  •  . 

m 

In  subcritical  case  is  not  necessary  to  put  coupling  conditions  but  in“ 
supercritical  case  is  absolutely  necessary  to  put  this  kind  of  conditions  in 
order  to  realize  the  continuity  of  the  velocity  field.  In  this  order  we  will  use 
Beavers-Joseph  postulate  according  to  which  the  slip  velocity  at  the 
permeable  interface  differs  from  the  mean  filtration  velocity  within  the 
porous  matrix,  and  the  shear  effects  are  transmitted  into  the  body  of  the 
material  through  a  boundary-layer  region.  Across  this  boundary  the  velocity 

changes  rapidly  from  this  value  at  the  interface  to  the  Darcy’s  law  value 
given  by 


K'  dp 


in  porous  medium  (like  in  Fig.3.). 


i 

L  Z 

U  \ 

. . -  ^ 

1  fluid  layer 

w 

_ V 

Um 

- ► 

porous  region 

w 

Fig.3.  Physical  model  for  Beavers-Joseph  slip  condition 


Beavers  and  Joseph  also  assumed  that  the  slip  velocity  for  the  free  fluid  is 
proportional  to  the  shear  rate  at  the  permeable  boundary  and  is  related  to  the 
slip  velocity  of  the  exterior  flow  by  an  adhoc  boundary  condition 


(16)  f\ 

dz 


=  cD.(u-u„) 


where  o+  is  the  boundary  limit  point  from  the  exterior  fluid,  and  O  is  a 
constant  which  depends  on  the  properties  of  the  fluid  and  the  porous  matrix, 

but  is  independent  of  spatial  coordinates.  They  showed  that  0  = 

where  a  is  a  dimensionless  quantity  depending  on  the  material  parameters 
which  characterizes  the  structure  of  the  permeable  material  within  the 
boundaiy  region.  Therefore  the  Beavers-Joseph  slip  condition  at  the 
interface  can  be  written  as: 

(17)  “  (u_u^) 

dz  Vic 

where  u  and  du/dz  are  evaluated  at  the  interface  o+,  while  Um  is  evaluated  at_ 
some  small  distance  d  given  by 

(18) 

a 


away  from  the  plane  z  =  0  in  porous  side.  Since  the  parameter  a  is  obtained 
from  a  dimensional  analysis  and  is  not  defined  mathematically  it  needs  to  be 
evaluated  from  the  experimental  data.  Beavers  and  Joseph  conducted 
experiments  to  determine  the  validity  of  the  above  model  and  reported  a  for 
three  kinds  of  foametal  and  two  kinds  of  aloxite  matrices  (6.5  x  10‘‘°  <  K'  < 
8.2  X  10'*  m^).  Their  data  also  demonstrated  that  a  is  almost  independent  of 
the  fluid  viscosity  p  (see  [48]).  Although  the  Beavers-Joseph  condition  was 
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first  presented  as  an  empirical  result,  a  certain  amount  of  theoretical 
foundations  was  later  provided  and  several  experimental  investigations 
supported  it.  In  any  case  Beavers-Joseph  condition  permits  to  find  the  small 
distance  d  where  we  make  the  coupling  between  the  slip  velocity  and  the 
mean  filtration  velocity  with  formula: 


Uc  = 


z  +  u  for  -  d  <  z  <  0 . 


In  the  paper  [12]  Jones  argued  that  the  Beavers-Joseph  condition  was 
essentially  a  relationship  involving  shear  stresses  rather  than  just  the  velocity 
shear,  and  suggested  a  generalized  form  (17)  as 


dii  dw  a  ,  X 


Jones  -  modified  Beavers  -  Joseph  condition 


where  v  =  (u,  v,  w). 


It  IS  understood  that  u,  — ,  —  are  evaluated  at  the  interface  z  =  0+,  while  Um 

(JZi  C/X 

is  evaluated  at  some  small  distance  d  given  by 


(20)  d  = 


Taylor  in  [49]  observed  that  the  Beavers-Joseph  condition  can  be  deduced  if 
the  Bnnkman  flow  model  is  used  in  a  region  very  close  to  the  fluid-porous 
layer  interface  while  the  flow  in  the  bulk  of  the  porous  medium  is  governed 
by  Darcy  law.  The  idea  was  extended  by  Neale  and  Nader  in  [50]  who 
showed  that  for  the  problem  of  flow  in  a  channel  bounded  by  a  thick  porous 
wall,  the  Brinkman  equation  yields  the  same  solution  of  that  obtained  from 

the  Darcy  equation  provided  a  is  taken  as  where  p'  is  the 

effective  viscosity  for  the  fluid/saturated  porous  medium. 

Ross  in  [51]  developed  an  equation  for  the  average  fluid  velocities  in  terms 
of  gradients  of  average  pressure  and  velocities  and,  indirectly,  the  geometiy 
of  the  porous  interface  by  considering  the  flow  of  inertialles  fluid  through 
the  pores.  Using  an  averaging  procedure,  he  obtained  the  interfacial  slip 
condition  for  the  Beavers-Joseph  flow  problem  as: 

(21)  u  =  + 

p  dx  dz  (12^ 

The  term  of  the  left  hand  side  and  the  first  two  terms  on  the  right-hand  side 
of  this  equation  are  identical  to  that  in  equation  (17)  except  that  the  first  term 
(equation  (21))  represents  the  Darcy  flow  at  the  interface  z=0  and  not  at  z<0 


as  in  the  Beavers-Joseph  model.  Therefore  neglecting  the  second  order 
viscous  term  in  (21)  and  considering  the  Darcy’s  flow  at  z<0  we  obtain  the 


Beavers-Joseph  model  with  a  = 


Ross  employed  a  scale  analysis,  to 


show  that  all  the  three  terms  in  equation  (21)  are  of  equal  order-of- 
magnitude,  and  the  second  order  viscous  term  must  be  retained  for  the  slip 
condition  to  be  valid.  A  generalized  slip  condition  for  an  anisotropic  porous 
medium  was  presented  as: 


(22) 


Ui=- 


5xj 


Ljpq-^  +  Kij-Njp-V^Up, 


where  Ky  =p-Kij  and  Ky  is  the  permeability  tensor  and  Ljpq  ,  Njp  are 
expressed  in  terms  of  line  and  area  integrals. 

Note  that  in  the  above  results,  the  porous  medium  is  dense  and  of  large 
thickness  as  compared  to  the  interface  boundary  layer. 

Rudraiah  in  [26]  examines  the  effect  of  porous  layer  thickness  on  slip 
velocity  at  the  interface,  particularly  when  the  porosity  is  large.  Modified 
Beavers-Joseph  expressions  were  obtained  for  two  cases:  a)  porous  layer 
bounded  below  by  a  thick  layer  of  static  fluid,  and  b)  porous  layer  bounded 
below  by  an  impermeable  rigid  wall.  For  case  b)  the  slip  condition  was 
found  as: 


sinh(5  •  h') 


+  (u-Un,)-coth(5h') 


where  h'  is  the  thickness  of  the  porous  layer,  5  =  and  V  is  a- 

viscosity  parameter.  Equation  (23)  reduces  to  the  Beavers-Joseph  condition 
as  h'->+oo  and  a  comparison  with  equation  (17)  shows  that  X'  =  a^,  which 
agrees  with  the  finding  of  Weale  and  Nader  in  [50]  if  X'  can  be  identified  as 
the  viscosity  ratio  p'/p.  It  is  understood  that  in  equation  (21),  u,  du/dz  are 
evaluated  at  the  interface  z  =  0+  while  u^  is  evaluated  at  some  small  distance 
d  given  by 


(24)  d  = 


Vic 


a 


-i-i 


coth(5  •  h')  +  ^  ^ 


u 


m 


sinh(5  •  h')  u  -  u , 
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Returning  to  the  mushy  zone,  this  is  located  in  the  region 
^  ~  t(^  i>Xj,X3)€R  |0<x,  +  <R  ,0<X3  <h'|  and  the  fluid  layer  in 
(he  region  n'=  {(x,,x,,x  Js  R=  |  0  <  xj  +  xj  <  R>,h'<  x,  <  L.}.  The 
permeable  interface  is  located  on  the  level  X3  =  h'.  Due  to  the  small  thickness 
and  large  porosity  of  the  mushy  zone  we  will  use  Rudraiah  slip  condition  for 
assuring  the  continuity  of  the  velocity  field. 

Once  global  velocity  field  v  and  pressure  field  p®  in  the  melt  are  determined 
from  the  two  phase  natural  convection  problem  described  by  equations  (10), 
(1 1),  (12)  in  super  critical  case  we  use  the  pressure  field  p°  in  order  to  find’ 
with  formula  (14),  the  Darcy’s  flow  <  v>.  With  this  flow  <  v>  and  the 
global  velocity  field  v  using  (22)  we  find  the  small  distance  d  where  we 
make  the  coupling  between  the  global  flow  v  and  Darcy’s  flow  <v>. 
Therefore  we  will  have  the  Darcy  flow  just  in  the  region  0  <  X3  <  d  <  h’. 


L3  The  transport  equation  of  the  dopant  in  mushy  zone  when 
convection  and  diffusion  balance  each  other 

Using  homogenization  method  in  the  case  when  convection  and  diffusion 
balance  each  other  (at  the  global  level)  we  can  obtain  (see  [52])  the 
following  transport  equation  for  a  solute  in  a  periodic  porous  media: 


6c 

- j_ 

a 


<Vi  > 


6c 

6x; 


i=l  j=l 


6^c 

6xj6xj 


In  this  equation  c  —  c(t,x)  is  the  solute  concentration,  <Vi>  =  <Vi>(xi,  X2,  X3), 
i  -  1 ,2,3  are  the  components  of  the  fluid  velocity  given  by  Darcy's  law  and" 
the  coefficients  Djj  are  the  components  of  the  diffusion  (dispersion)  tensor. 
These  components  are  given  by  the  formula 


(26)  Dii=D 


5„+  '  fedy 


=  D 


5jj  + 


In  the  above  formulas  D  is  the  diffusion  (dispersion)  coefficient  of  the  solute 
and  Xi  is  the  solution  of  the  following  boundary  value  problem 
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(27) 


AyXi  =0  inYp 

n-VyXi=-ni  onT 
<Xi  is  Y  periodic 


<Xi>=|Y| 


Jxidy 

Yf 


=  0 


for  i  =  1,  2,  3. 


Multiplying  equation  AyXj  =0  by  Xj  and  integrating  on  Yp  we  obtain  for 
the  components  Dy  the  following  formula 


(28)  Di:=D, 


Su-j^fvXi-VXjdy 


and  consequently  Dy  is  symmetric  (Dy  =  Dj;). 
In  addition  we  have 


(29)  Dij=D.<VBi,VBj>=D.jl 

where  Bj  =  -Xj  -  Yj  and  consequently  the  tensor  Dy  is  positive. 


If  Ys  is  symmetric  to  any  coordinate  planes  then  Dy  =  0  for  i  j  and 
consequently  the  dispersion  tensor  is  diagonal. 

When  the  solid  inclusions  are  small  with  respect  to  the  distance  between  two 
neighbouring  inclusions  using  results  from  [46]  and  [47]  in  [53]  was  proved 
that  for  symmetric  solid  inclusions  we  have 


(30)  Dii=Di 


where  (p  is  the  volume  concentration  of  solid  inclusions  and  is  assumed  to  be 
small.  The  functions  Xi  are  solutions  of  the  external  Neumann  problem: 

(31)  AXi=0  inR^Ys 


(32)  AXj*n=-nj  onT 
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In  [54]  we  have  find  the  transversal  (Di  1=022)  and  longitudinal  (D33) 
diffusion  coefficients  for  the  prolate  spheroids.  They  are  plotted  as  function 
of  the  excentricity  on  the  following  figures. 


Fig.4. 

Therefore  the  dopant  transport  in  mushy  zone  is  given  by  (25)  where  Dy  are 
given  by  (30)  and  <  Vj  >  are  the  component  of  the  velocity  field  of  the 
global  convection  flow,  Brinkman’s  flow  or  Darcy  flow. 

Steady  state  solutions  of  equation  (25)  are  given  by  equation 


-i 

Y  i=l 


<V:  > 


i=l  j=l  C7X,C7Xj 


5x 


This  equation  generalizes  the  steady  state  mass  balance  equation 

(34)  — (v-Vc)=V^c 
Pr 

used  in  [20]  for  the  determination  of  the  dopant  field  in  the  melt.  In  equation 

(34)  the  Schmidt  number  Sc  =  p/3)  and  the  Prandtl  number  Pr  =  p/ttL-  3)  is 

the  diffusivity  of  Ga  in  Ge,  aL  is  the  thermal  diffusivity  in  liquid.  Boundary 
conditions  used  in  [20]  for  solving  (34)  are 

(35)  N-Vc  =  J^(e,-N).(l-k)c;  0<r<A.  z  =  h(r) 

(36)  —  =  _^(c~l)  0<r<A,  z  =  0 

fe  Pr 

(37)  —  =  0;r  =  0,A,  0<z<l 
or 
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where  z  =  h(r)  represents  the  location  of  the  melt/solid  interface,  A  =  R/L,  R 
is  the  radius  and  L  is  the  length  of  ampule,  N  is  the  unit  vector  normal  to 
the  interface,  Pe  =  Vs-L/at  is  the  Peclet  number,  k  is  the  segregation 
coefficient  of  Ga  in  Ge. 

Equations  (35),  (36)  express  conservation  of  mass  at  the  melt/solid  interface, 
and  the  fictitious  "inlet"  at  the  melt  end  (z  =  0)  of  the  ampule  respectively. 
Equation  (37)  is  the  no  flux  condition  valid  at  the  centerline  and  sidewalls  of 
the  ampule. 

Once  the  velocity  field  in  the  melt  and  the  shape  of  melt/solid  interface  were 
determined  from  the  two  phase  natural  convection  problem:  equations  (34)  - 

(37)  are  reduced  to  a  linear  set  to  be  solved  for  concentration  distribution 
through  the  melt. 

Returning  to  equation  (33)  boundary  conditions  for  solving  this  equation  can 
be  formulated  starting  from  the  location  of  the  mushy  zone. 

We  have  already  supposed  that  the  mushy  zone  is  located  in  the  region 

(38)  n  =  |(x,,X2,X3)|xf  +X2  <R^  and  0<X3  <  h'j 


and  from  (35)  it  follows  that  we  have  equation 

(39)  osr<R, 


dx 


Pr 


=  0 


which  expresses  conservation  of  mass  at  the  melt/solid  interface. 

From  equation  (37)  which  expresses  the  no  flux  condition  at  the  centralline 
and  sidewalls  it  follows 

(40)  —  =  0,  r  =  0,R;  0<X3<d 
or 


Equation  (36)  must  be  substituted  by  a  coupling  condition  at  the  level  X3  =  d 
between  the  concentration  field  c  defined  by  (34)  -  (37)  and  the" 
concentration  field  c  which  satisfies  (33),  (39)  and  (40)  for  X3  <  d.  We  can 
realize  this  coupling  putting 

(41)  c(xi,X2,d)=c(x,,X2,d) 

Once  the  velocity  field  in  "mushy  zone"  are  determined  (from  equations 
(10)  -  (12)  or  from  (13)  or  from  (14))  and  c  determined  from  (34)  -  (37)  we 
can  solve  (33),  (39),  (40),  (41)  and  we  have  the  concentration  distribution  in 
mushy  zone.  This  dopant  field  adjacent  to  the  melt/solid  interface  is  a  new 
more  rafined  approach  of  the  real  dopant  field  in  this  region. 
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1.4  The  transport  equation  of  the  dopant  in  mushy  zone  when 
convection  dominates  diffusion 


Using  homogenization  method,  in  the  case  when  convection  dominates 
diffusion  at  global  level,  we  can  obtain  (see  [30]),  [31]  and  [55])  the 
following  transport  equation  for  a  solute  in  a  periodic  porous  media  Q. 


In  the  above  equation  c  —  c(t,x)  is  the  solute  concentration, 

<Vi>  =  <Vi>(xi,  X2,  X3),  i  =  1,  2,  3  are  the  components  of  the  velocity  field 
and  the  coefficients  Dy  are  the  components  of  the  diffusion  (dispersion) 
tensor.  These  components  are  given  by  the  formula 


-raft" 


In  the  above  formula  D  is  the  diffusion  coefficient  of  the  solute  and 
Xi  ~  Xi  (^5  y)  is  the  solution  of  the  following  boundary  value  cell  problem 
AyXi -V- VyXi  =Vi-<Vj  >  inYp 
■  ^yXi  =  “Hj  on  r 

(44)  .  Xi  is  Y  periodic 

j3Cidy  =  0 

V  is  the  velocity  field  defined  by 

(45)  v(x,y)=-^^.vJ 

j=i 

where  are  defined  by  (8),  p®  is  the  pressure  field  in  Q  and  <Vi>  =  <Vi>(x) 
is  given  by 

(46)  <Vi  >(x)  =  jl  Jvi(x,y)dy. 

For  a  constant  gradient  pressure  field  Dy  are  constants  and  due  to  the  identity 

M  j=i  c^i-c^Xj  i^j  j^j2  axi-axj 
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in  equation  (42)  we  can  use  only  the  symmetric  part  D-  = 


the  diffusion  (dispersion)  tensor.  This  symmetric  part  is  given  by  one  of  the 
following  three  equivalent  formulas: 


(48)  D)j=D 


|(xiVj  +XjVi)tiy 


Sij  +XjVi)dy 

3yk 

where  Bi  is  defined  by 
(51)  B,=-Xi-yi. 

If  the  mean  flow  is  directed  along  the  X3  axis  then  Dy  is  diagonal  with  two 
independent  components  that  are  the  longitudinal  Dl  and  transverse 
dispersivities  Dt;  Dl  =  D33  and  Dj  =  Du  =  D22  (see  [17],  [18]  and  [26]).  For 
any  other  flow  direction  in  the  plane  Oyiy2,  there  are  four  nonzero 
independent  dispersivity  coefficients  Du,  D22,  D33  and  D12  =  D21;  D13  =  D31 
=  D23  =  D32  =  0  (see  [17],  [18]  and  [26]). 

To  calculate  the  dispersion  tensor,  one  must  solve  the  canonical  cell  problem 
defined  by  (44).  Lee  Sun  and  Mei  in  [18]  have  shown  that  the  cell  boundary 
value  problem  for  convective  diffusion  of  heat  can  be  recast  as  a  variational 
principle,  the  case  for  solute  transport  being  a  special  case.  Computed  values 
of  longitudinal  and  transverse  dispersivity  coefficients  Dl  and  Dj,  for  a 
cubic  array  of  Wigner-Seitz  grains,  for  Peclet  number  up  to  300  for  Dl  and 
200  for  Dt,  and  for  two  porosieties  7tY"=  0.38  and  0.5  are  shown  in  the" 
following  figures. 


(49)  Di=D 

(50)  D(:=D 
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To  conform  with  experimental  literature  the  abscissas  are  the  Peclet  numbers 
defined  in  terms  of  the  mean  flow  velocity  averaged  over  Yp  only,  i.e. 
<V|>-?i/71y-D  =  Pe/Tiy.  In  Fig.5(a),  the  longitudinal  dispersivity  is  also 
compared  with  the  measured  data  [57]  for  simple  cubic  packing  of  uniform 
spheres  and  the  calculations  [58]  for  a  simple  cubic  lattice  of  uniform 
spheres  with  Tty  =  0.48,  0.74  and  0.82.  The  results  for  Tty  =  0.48  by  Koch  et 
al.  [59]  based  on  an  approximate  analysis  for  dilute  concentrations  are  also 
included.  All  are  in  qualitative  agreement  for  Dl.  From  the  numerical  results 
for  n  —  0.38  and  0.5  we  see  that  for  small  Peclet  numbers  where  molecular 
diffusion  is  dominant  in  microcell,  the  diffusivity  is  greater  for  the  larger 
porosity.  The  reason  is  that  the  cross-sectional  area  through  which  a  passive 
solute  can  diffuse  increases  with  porosity.  It  is  always  less  than  unity  in  the 
diffusion-dominated  regime  (at  microlevel)  since  the  presence  of  solid  grains 
reduces  the  diffusive  flux  of  solute.  For  large  Pe  numerical  results 
concerning  Dl  for  Wigner-Seitz  cell  as  well  as  those  by  Salles  et  al.  for- 
uniform  spheres,  are  consistent  with  the  experimental  measurements  [58] 
and  the  analytical  theory  for  dilute  spheres  [59],  all  showing  that  Dl 
increases  with  ^e) ,  when  the  mean  flow  is  parallel  to  a  lattice  axis.  Recall 
from  [59]  that  if  the  flow  is  inclined  to  a  lattice  axis,  Dl  may  vary  linearly 
with  Pe.  In  contrast  to  the  case  of  small  Pe,  the  dependence  on  porosity  is 
now  reversed,  and  the  dispersivity  increases  with  decreasing  porosity. 
Heuristically  this  is  because  the  velocity  gradient  in  the  pores  increases  as 
porosity  decreases  and  therefore  enhances  microscale  mixing. 

The  transverse  dispersivity  Dy  computed  for  Wigner-Seitz  grains  is  plotted 
m  Fig.5(b)  for  Pe  <  300.  The  qualitative  trend  is  the  same  as  Dl  except  that 


22 


it  is  less  than  Dl  by  roughly  two  orders  of  magnitude.  Mauri  [60]  also  finds 
analytically  for  a  dilute  lattice  of  uniform  spheres  that  Dj  is  proportional  to 
(Pe)^,  for  small  Pe  and  is  eight  times  smaller  than  Dl.  In  contrast  Koch  et  al. 
in  [59]  predict  that  Dj  remains  almost  constant  in  Pe  for  very  large  Pe.  There 
are  no  reliable  measurements  for  Dt  for  a  regular  array  of  spheres.  Some 
experimental  data  on  Dt  for  natural  granular  media  are  shown  in  the 
following  figures  (see  paper  [61]  -  [65]).  Although  scattered,  each  individual 
data  set  exhibits  linear  dependence  on  Pe  as  Dl. 


Pt/n 


£  Cumpanson  of  lar  i;fuin  wuh  dau  for  naiufsl  grarudar  media. 


Steady  state  solutions  of  equation  (42)  for  a  constant  gradient  pressure  field 
are  given  by 

(-)  =i2:o. 


i=i  j=i 


0Xj5xj 
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where  D-j  is  given  by  one  of  the  three  equivalent  formulas  (48),  (49),  (50) 
and  <  Vj  >  are  the  components  of  the  velocity  field. 

Equation  (52)  generalizes  the  mass  balance  equation  (34)  used  in  [20]  for  the 
determination  of  the  dopant  field. 

For  solving  equation  (52)  we  can  use  the  boundary  conditions  (39)  (40)  and 
(41). 

Once  velocity  field  in  mushy  zone  is  chosen  we  can  solve  the  boundary 
value  problem  (52),  (39),  (40)  and  (41)  finding  the  dopant  distribution  in 
mushy  zone.  This  dopant  distribution  adjacent  to  the  melt/solid  interface  is  a 
more  rafined  approach  of  the  real  dopant  field  in  this  region. 


1.5  Stability  and  segregation 


Solving  boundary  value  problems  (33),  (39),  (40),  (41)  or  (52),  (39),  (40) 
and  (41)  we  obtain  dopant  distributions  which  are  maybe  in  a  better 
agreement  with  experimental  data  and  reflect  the  real  radial  segregation.  But 
the  main  question,  which  rests  to  be  solved,  is:  in  which  way  we  can 
obtain  in  a  steady  state  regime  a  semiconductor  crystal  without  radial 
segregation? 


In  our  opinion  in  order  to  obtain  in  a  steady  state  regime  a  semiconductor 
crystal  without  radial  segregation,  the  dopant  concentration  in  the  mushy 
zone  most  be  independent  on  t,  x,,  X2.  In  this  sense  see  also  [66]  concerning 
experiments  conducted  on  board  LDLG  vehicles:  "the  dopant  field  close  to 
the  growing  surface  must  be  uniform  prior  to  low-g  growth".  A  dependence 
on  xi,  X2  may  lead  to  segregation.  This  means  that  in  mushy  zone  the 
concentration  c  of  the  dopant  may  depend  only  on  X3,  c  =  c(x3). 

Using  this  kind  of  functions  in  the  dopant  transport  equation  we  find  that  c 
satisfies  the  equation. 

(53)  tty 


dx 


y=<V3  > 


dX' 


Now  if  <V3>  depends  on  Xj  and  X2  then  equation  (38)  has  only  constant 
solutions.  This  means  that  in  the  mushy  zone  we  must  have  a  constant 
concentration  of  the  dopant.  This  corresponds  to  the  Griffin  and  Motakef 
recommendation  made  in  [66]  but  does  not  satisfy  the  boundary  condition 
(39).  If  <V3>  depends  only  on  X3  we  have  several  solutions  for  the  equation 
(53)  and  to  each  of  these  corresponds  a  dopant  distribution  in  mushy  zone. 
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without  radial  segregation.  Obviously  we  must  choose  solution  which 
satisfies  boundary  condition  (39),  (40)  and  (41).  ^ 

In  particular  if  in  mushy  zone  we  have  Ky  =  0  for  i  j  and  =  0  then 

<V3>  =  0  and  we  have 

(54)  Co(x3)=ki +k2X3 

where  ki,  k2  are  two  constants  to  be  find  from  (39)  and  (41).  The  boundary 
condition  (40)  is  satisfied  for  any  constants  kj  and  k2. 

If  in  a  steady  state  regime  the  dopant  concentration  in  mushy  zone  is  given 
by  a  solution  of  ecjuation  (53)  then  this  rests  the  same  all  the  time,  till  the 
regime  does  not  change  and  we  will  obtain  a  semiconductor  crystal  without 
radial  segregation. 

Now  the  main  question  is:  starting  from  a  certain  steady  state  regime 
and  from  the  corresponding  dopant  distribution  in  which  “way”  we 
succeed  to  realize  a  new  steady  state  for  which  the  dopant  distribution  is 
given  by  a  solution  of  equation  (53)? 

In  order  to  give  a  partial  answer  to  this  question  we  will  try  to  explore  the 
region  of  attraction  (region  of  stability)  of  a  solution  of  equation  (53).  This 
means  that  for  a  solution  Cq  “  Co(x3)  of  equation  (53)  we  search  perturbations 
which  tend  to  0  for  t  tending  to  +oo.  If,  for  a  solution  cq,  the  class  of  this  kind 
of  perturbations  is  sufficiently  large  then  we  can  try  to  find  a  perturbation 
such  that  the  sum  of  Cq  ~  Co(x3)  and  the  perturbation  at  t=0  coincides  with  the 
dopant  distribution  corresponding  to  the  starting  steady  state.  If  we  succeed 
to  find  such  a  perturbation  then  we  have  find  the  “way”. 

For  the  beginning,  we  consider  perturbation  of  the  form 

(55)  Cj(xi,t)=c,(xi)e'^ 

and  we  replace  c  in  the  transport  equation  (25)  by 


(56)  c  =  Co(x3)+Ci(xi)-e‘^. 


We  find  that  Cj  =  Ci(xi)  satisfies  the  equation 


(57)  a-c, +  — <vi  >-^  =  D 


dci  _  d^Cj 


dxi 


j  2 

dxi 
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If  <Vi>  depends  on  X2  and  X3  then  c,  =  0  is  the  only  solution  of  equation 

(57) .  But  Cl  =  0  is  not  a  real  perturbation.  If  <vi>  depends  only  on  Xi  then 
equation  (57)  has  several  solutions.  In  particular,  for  any  real  a  and  8  >  0 
equation  (57)  has  solutions  with  the  following  property: 

(58)  sup|c,(x,)|<E 

where  the  supremum  is  considered  on  the  domain  0  <  Xi  <  R.  This  means 
that  the  above  considered  solution  Cq  =  Co(x3)  is  not  stable  with  respect  to  any 
perturbation  of  the  form  (55). 

If  we  impose  to  the  perturbation  (55)  to  satisfy  the  boundary  condition  (40) 
then  we  obtain 

(58)  —^  =  0  forxj  =0andxi  =R. 

QX  j 

For  the  Sturm-Liouvill  problem  (57),  (58)  there  exist  a„  and  c,„  such  that 
0  >  a,  >  a2  >  G3  >  ...  >  -00,  c,n  =  c,n(xi)  0  and  Gn,  c,n  satisfy  (57),  (58).  This 
means  that  the  concentration  profile  Co  defined  by  (53)  is  asymptotically 
stable  with  respect  to  the  perturbation 

(59)  Ci„(x,,t)  =  c,„(x,>'^"'. 

Now  if  we  study  the  stability  of  the  same  concentration  profile  with  respect 
to  a  perturbation  of  the  type 

(60)  C2(x2,t)  =  C2(x2)-e°' 

we  can  establish  similar  facts. 

If  we  investigate  the  stability  of  the  concentration  profile  with  respect  to  a 
perturbation  of  the  type 

(61)  Ci2(xj,X2,t)  =  Ci2(xi,X2)-e'^ 

we  will  find  other  perturbations  with  respect  to  which  the  dopant 
concentration  profile  cq  =  Co(x3)  is  asymptotically  stable. 

This  stability  results  prove  that  there  are  time  dependent  regimes  for  which 
the  dopant  distribution  in  mushy  zone  tends  to  a  steady  dopant  distribution 
profile  without  radial  segregation.  This  kind  of  facts  are  encouraging  but 
don't  mean  great  things  until  we  have  not  find  the  whole  family  of 
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perturbations  with  respect  to  which  the  dopant  concentration  profile  Cq  = 
Co(x3)  is  asymptotically  stable. 

1.6  The  optimal  boundary  control  for  a  solution  without  radial 
segregation 


Let  us  consider  in  the  flow  domain  M  (the  melt)  the  natural 
convection  governed  by  the  equations: 

(62)  -  Pr-  V  +  (vV)v  +  Vp  =-  Ra  •  Pr-  0  •  es 

(63)  -V^0  +  v-V0  =  O 

(64)  V  •  V  =  0 

and  the  dopant  dispersion  governed  by  the  equation 

(65)  -Pr-V^c  +  Sc-v- Vc  =  0. 


We  assume  that  the  boundary  conditions  are  the  followings: 

(66)  V  =  0  on  5M 


(67)  ^  =  a(0-T) 
on 


on5M 


(t  is  the  temperature  of  surrounding  medium) 
dc 

(68)  — =-  =  y*c-n3  on5M 

dn 


ns  =n-e3. 


The  aim  is  to  find  the.  boundary  controls  x,  which  give  us  a  desired 
field  of  concentration  Cd  uniform  in  neighborhood  of  the  melt/solid  interface. 
This  means  to  minimize  the  cost  functional  defined  by 

(69)  J(c,x)=^  J]c-Cd|  dx  . 

M 

We  will  present  some  results  which  are  related  and  can  be  used  in 
order  to  solve  the  above  boundary  control  problem. 

Abergel  and  Casas  in  the  paper  [67]  consider  a  more  general  system  as  (62)- 
(64)  with  boundary  conditions  similar  to  conditions  (66)-(67)  and  solve 
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boundary  control  problems  which  consist  in  minimizing  turbulence  caused 

by  the  heat  flux  on  the  boundary  i.e.  to  minimize  the  cost  functional 
involving  the  turbulence,  defined  by 

(70)  j(v,T)=i|vxq'dx  +  ||9-f ds. 

M  dM 

Capatina  and  Stavre  in  [69]  consider  also  a  more  general  system  as  the 
system  (62)-(64)  with  the  boundary  conditions  (66)-(67)  and  solve  the 
boundary  control  problems  which  consist  to  find  the  control  x  which  give  as 
a  desired  field  of  temperature  Gj ;  i.e.  to  minimize  the  cost  functional: 

(71)  J(e,t)=ye-e,fdx  . 

M 

[70]  Capatina  and  Savre  give  a  numerical  treatment  of  the  above 
boundary  control  problem.  Finite  element  approximation  of  the  optimality 
system  is  defined.  The  convergence  of  the  proposed  algorithms  for  solving 
the  discrete  problem  is  proved.  The  analysis  of  the  numerical  results  and 
their  physical  meaning  are  discused. 

Finally  Capatina  and  Stavre  in  [71]  solve  an  optimal  control  problem  in 
biconvective  flow;  “biconvection”  being  convection  caused  by  a 
concentration  gradient,  not  by  a  temperature  gradient.  They  find  the  mean 
values  a  of  concentrations  which  lead  to  a  given  concentration  field  Ch  They 
consider  the  cost  functional 

(72)  J(a,c)  =  ^  J]c-Cd|  dx  +  ~a^ 

M 

the  optiiriBl  control  probicni  ss  follows* 

(73)  min{j(a,c)|(a,c)eT} 

where  T  is  the  nonempty,  weakly  closed  set: 

(74)  T  =  {(a,  c)  6  [0, 00)  X  H  ‘  (q)|  3  v  such  that  (v,  cjsatisfies  (62)  -  (64) j 

The  physical  relevant  term  in  (72)  is  j  |c-c,|'dx  which  provides  an 

M 

estimate  of  the  difference  between  the  component  c  of  (v,c),  and  a  given 
configuration  Cj  of  concentration. 

This  result  is  very  interesting  in  the  case  of  CdZnTe  crystal  growth  for 
which  the  solute  Rayleigh  number  is  Ras=2.02xl0*  (for  GeGa  the  solute 
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Rayleigh  number  is  Ras=0)  and  consequently  equation  (62)  must  be  replaced 
by  the  equation: 

(75)  -  Pr- v  + (vV)v  +  Vp  =- (Ra  •  Pr- 0  +  Ra^  •  Pr- c)e3 


7.  7  Summary  and  Recommendation 


1.  By  a  synthesis  of  the  results  concerning  the  mass  transport  in  porous 
media  we  have  identified  two  transport  equations  (25)  and  (42)  which 
generalize  the  connective-diffusive  mass  transport  equation  (4)  used  in  [20] 
and  in  general  for  the  determination  of  the  dopant  field. 

2.  Locally,  in  the  precrystallization  zone,  which  can  be  physically  localized, 
equations  (25)  and  (42)  describe  with  more  accuracy  the  dispersion  of  the 
dopant  due  to  the  structure  which  already  exists  in  this  region  and  the 
interaction  between  complicated  global  flow  patterns  and  the  convective- 
diffusive  mechanism  which  exists  simultaneously  in  this  region.  Equation 
(25)  describes  the  dispersion  process  when  convection  and  diffusion  balance 
each  other,  equation  (42)  describes  the  process  when  convection  dominates 
diffusion. 

3.  Solving  numerically  boundary  value  problems  (35),  (39),  (40)  and  (41) 
respectively  (52),  (39),  (40)  and  (41)  for  Ga  in  Ge  we  hope  to  obtain  better 
agreements,  like  Chang  and  Brown  in  [20],  with  experimental  data,  in  the 
case  of  steady  state  regime. 

4.  For  this  reason  we  recommend  the  implementation  of  this  transport 
equation  in  M.A.S.T.R.A.P.  computational  model. 

5.  It  must  be  noted  that  equations  (25)  and  (42)  do  not  take  into  account  the 
interaction  of  the  dopant  field  with  the  crystal.  A  kind  of  interaction  appears- 
just  in  boundary  condition  (39).  In  this  context  we  can  think  to  improve 
equations  (25)  and  (42)  by  the  addition  of  a  new  term  which  expresses  this 
interaction.  Something  like  the  terms  which  express  absorption  and  reaction 
in  equation  of  Homung  and  Yager  in  [52]. 

6.  The  stability  calculus  made  for  the  steady  concentration  profiles  without 
radial  segregation  suggests  that  rotating  the  ampule  around  the  longitudinal 
axis  we  can  reduce  <V3>  and  obtain  a  quasi  Couotte  flow  in  mushy  zone, 
reducing  in  this  way  the  radial  segregation. 

7.  The  optimal  boundary  control  calculus  is  just  a  beginning  and  suggest  in 
how  way  we  can  reduced  the  radial  segregation. 
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2.  On  the  determination  of  the  concentration  of  the  solute  in 
the  case  of  crystal  growth  by  hydrothermal  method 


2.1  Introduction 

We  consider  the  prototype  vertical  hydrothermal  system  for  obtaining  quartz 
single  crystals  presented  in  Fig.  1 . 


tuld  zom 


/:onc 


hot  -/jono 


ooooooooo 


a)  autoclave 


b)  temperature  distribution 


Fig.l.  Geometry  of  the  prototype  vertical  hydrothermal  system 

Using  the  transport  model  of  solutes  in  a  porous  medium  participating 
in  a  dissolution  reaction  in  general  not  in  equilibrium;  developed  by  Knaber, 
Duijn  and  Hengst  in  [72];  for  the  lower  region,  in  the  porous  bed,  if  the 
charge  distribution  is  constant,  we  find  travelling  waves.  The  travelling 
wave  in  fact  exhibits  a  sharp  dissolution  front  and  is  obtained  in  a  nearly 
explicit  manner  for  two  spatially  one-dimensional  flow  regimes  with 
constant  water  content,  bulk  density,  pore  velocity  q  and  diffusion 
coefficient  D.  Also  for  the  limit  cases  of  equilibrium  reaction  or  no 
dispersion,  travelling  waves  are  obtained  under  the  same  conditions,  but 
with  different  qualitative  properties  (see  also  [73]).  It  must  be  noted  that  the 
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spatial  domain  is  assumed  to  be  large  due  to  that  solutions  describe  the 
behavior  for  large  flow  domain. 

Chen,  Prasad,  Chatterjee  and  Larkin  in  a  recent  paper  [74]  since  the 
dissolution  and  growth  process  is  very  slow  neglect  the  effect  of  dissolution 
and  growth  and  consider  that  in  the  lower  region  (in  porous  bed)  the.  flow  is 
a  Darcy-Brinkmann-Forchheimer  flow  model.  Their  study  focuses  on  the 
flow  and  temperature  fields  without  nutrient  transport.  They  use  the 
following  nondimensionalized  governing  equations  in  the  lower  region 
(porous  region): 


(1)  Vu  =  S. 


(2)  — +  i(uV)-  =  -Vp  +  Gr-e-e, 
s  5t  s  s 

(3)  Se-  — +u-V0  =  — V(R,V0) 

^  ^  a  Pr  ^  ^ 


( 


+  VAVu - 


1 


■  + 


Fs 


Da  Da' 


u 


u 


where;  u  is  the  dimensionless  velocity;  Sm  the  mass  source  term;  8  the 
porosity;  p  the  dimensionless  pressure,  Gr  =  g  P-R  -(Th  -  Tc)/  Vf  is  the 
Grashof  number;  g  the  acceleration  due  to  the  gravity;  P  isobaric  expansion 
coefficient;  R  radius  of  cylindrical  autoclave;  Th  and  Tc  hot  and  cold 
temperature;  respectively  Vf  kinematics  viscosity;  0  =  (T  -  Tc)/(Th  -  Tc) 
dimensionless  temperature;  A  =  Peff^P  viscosity  ratio;  Da  =  K/R^  Darcy 
number;  K  permeability  of  porous  matrix;  Fs  =  b/R  Forchheimer  number;  b 
Forchheimer  coefficient;  Se  =  (Cp)efi/(Cp)f  specific  heat  ratio;  Pr  =  (v/a)f 
Prandtl  number;  Rk  =  kefi/k  ratio  of  thermal  conductivity;  a  thermal 
diffiisivity. 

The  flow  and  temperature  dynamic  governed  by  these  equations  correspontT 
to  an  isotropic  porous  media.  They  develop  a  three-dimensional  algorithm 
based  on  the  curvilinear  finite  volume  technique  and  a  non-staggered  grid 
layout  to  simulate  the  flow  and  heat  transfer  in  a  typical  autoclave  system. 
At  low  Grashof  numbers  an  axisymmetric  flow  pattern  and  at  high  Grashof 
numbers  a  three-dimensional  flow  pattern  is  predicted.  The  study  is  also 
extended  to  study  the  outset  of  oscillatory  flow  with  a  variation  in  the  porous 
bed  height.  In  the  upper  region  (cold  region)  Chen,  Prasad,  Chatter]  e  and 
Larkin  use  a  system  which  describe  the  natural  convection  without  the 
presence  of  the  seed,  governed  by  the  equations: 
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(4)  —  +  (vV)v  = -Vp  +  Pr  V^v  -  Ra  •  Pr- 0  •  e, 

ot  ' 

(5)  ^  +  v.v0  =  v-e 

(6)  Vv  =  0 

They  use  the  same  three-dimensional  algorithm  based  on  the  curvilinear 
finite  volume  technique  and  a  non-staggered  grid  layout  to  simulate  the  flow 
and  heat  transfer  in  the  cold  region. 

These  results,  for  the  first  time,  depict  the  possible  flow  patterns  in 
hydrothermal  system  that  can  have  for  reaching  consequences  on  the  growth 
process  and  crystal  quality.  ^ 

In  order  to  see  these  consequences  is  necessary  to  consider  the 
nutrient  transport  equation  in  homogeneous  media: 


/_\  dc  Sc  _  „ 

(7)  +_.v.Vc  =  V-c 

at  Pr 

and  to  transpose  it  in  a  porous  media  for  the  hot  zone  (lower  region  of  the 
autoclave)  and  to  analyze  the  accuracy  of  equation  (7)  in  the  neighborhood 
of  the  seed.  This  analysis  is  also  important  because  the  repartitioning  of  the 
nutrient  in  this  zone  influence  the  quality  of  the  grown  crystal.  In  this 
context  the  ‘diffusion-layer”  model  introduced  by  Noyes  and  Whitney  based 
on  a  simplified  form  of  equation  (7)  in  the  neighborhood  of  the  seed,  which 
importance  in  crystal  gro\^^h  from  solutions  has  been  stressed  by  Nernst  and 
which  was  used  by  several  authors  for  “flat”  or  “rough”  interfaces,  it  seems 
to  be  a  “gross  oversimplification”  of  the  interaction  between  complicated 
flow  patterns  and  solute  field.  Only  detailed  calculations  of  the  exact 
interaction  of  fluid  flow  and  solute  profile  adjacent  to  the  interface  can  be 
used  to  estimate  the  segregation  in  crystal  (see  [20]). 

Starting  from  this  also  for  the  vertical  hydrothermal  growth  system  we 
propose  the  substitution  of  the  “stagnant  film”  or  “diffusion  layer”  by 
the  “precrystallization  zone”  or  “mushy  zone”. 

This  zone  is  a  thin  bed  of  periodically  distributed  solid  inclusions,  which 
are  small  with  respect  to  the  distance  between  two  neighboring 
inclusions.  The  diffusion-layer  thickness  is  substituted  by  the  mushy 
zone  thickness,  which  can  be  determined.  In  this  way  the  diffusion-layer 
or  stagnant  film  masking  the  growing  crystal  is  substituted  by  a  thin 
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periodic  porous  medium  bounded  below  by  an  impermeable  rigid  wall 
which  is  the  fluid/solid  interface. 

If  the  size  of  solid  inclusions  is  smaller  like  a  critical  size  (see  [23])  the 
flow  in  “mushy-zone”  is  the  natural  convection  given  by  (4)-(6). 

If  the  size  of  solid  inclusions  is  critical  (see  [23])  then  the  flow  in  the 
mushy  zone  is  a  Brinkman  flow,  determined  by  the  global  pressure  and 
respecting  at  the  “mushy-zone”/”pure-fluid”  interface  the  Beaves- 
Joseph  or  Jones  modified  Beavers-Joseph  or  Rudriah  conditions  (see 
(241,  [25],  [26]). 

If  the  size  of  solid  inclusions  is  larger  like  the  critical  size  (see  [23])  the 
flow  in  “mushy-zone”  is  a  Darcy  flow,  determined  by  the  global 
pressure  and  respecting  at  the  “mushy-zone”/”pure-fluids”  interface 
Rudriah  conditions  (see  [26]). 

The  dispersion  mechanism  of  the  nutrient  in  mushy  zone  is  due  to  the 
diffusion  and  convection,  which  exists,  simultaneously  in  this  periodic 
porous  medium.  The  dispersion  mechanism  takes  into  account  the 
existence  of  solid  inclusions. 

This  dispersion  mechanism  of  the  nutrient  in  the  neighborhood  of  the 
interface  is  more  refined  as  the  dispersion  mechanism  in  diffusion-layer 
model  as  well  as  the  dispersion  mechanism  described  by  equation  (7) 
and  is  governed  by  a  new  convective-diffusive  equation  which  generalize 
equation  (7).  In  the  new  convective-diffusive  equation  we  use  the 
velocity  field  defined  by  equation  (4)-(6)  or  the  Darcy  flow  determined 
by  this  velocity  field  and  a  certain  periodic  micro-structure  similar  with 
the  micro-structure  of  the  crystal.  We  find  stationary  solutions  without 
segregation  for  the  new  convective-diffusive  equation  and  perturbation 
with  respect  to  this  solutions  are  assymtotically  stables.  We  establish 
also  the  controlability  of  the  stationary  solutions  without  segregation. 


2.2.  Porous  media  based  nutrient  transport  in  hot  zone 


In  the  following  we  will  describe  an  algorithm  for  the  determination  of 
the  transport  equation  of  the  nutrient  in  the  case  when  the  porous  medium  in 
the  hot  zone  is  periodic,  the  solid  obstacle  is  symmetric  and  convection- 
diffusion  balance  each  other; 
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1.  Let  be  L  (L  -  Hl,  Fig.2.  pg.23.  paper  Chen  et  al.)  the  size  of  the  porous 

medium  which  is  divisible  into  periodic  cubes  Q  of  dimension  £ 
presented  on  the  Fig.2. 


Fig.2. 

©  IS  the  solid  part;  Q'  =  Q  -  co  is  the  fluid  part;  dco  is  the  boundary  of®-  dQ 
is  the  boundary  of  Q.  ’ 

The  solid  part  ®  is  assumed  to  be  symmetric  with  respect  to  each  coordinate 
plane.  It  is  also  assumed  that  s  =  g  /L  «  1 . 

2.  In  dimensionless  variables  y  =  x/Ethis  cube  is  the  unite  cube  Y  presented 
in  Fig.3. 


Fig.3. 
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The  fluid  domain  Yp  and  the  solid  one  Ys  have  a  common  smooth 
boundary  F. 

3.  In  the  unite  cube  Y  we  have  to  solve  the  following  3-D  boundary 
value  problems 

AZi  =0  inYp 

(8)  VXi-n  =  -ni  onF  i  =  1,2,3 


Xi  is  Y  periodic 


Jxidy=o 

Yp 

where:  h  is  the  external  unit  normal  to  the  solid  boundary  F; 
^  =  the  condition  that  Xi  is  Y  periodic  means  that  we  have: 

Xi(-l/2,y2,y3)=Xi(l/2,y2,y3)  XjG' i -l/2,y 3)=  XjG' i.I/2,y 3) 
Xi6'l,y2 -1/2)=  Xjly^y  2,1/2)  i  =  1,2,3. 

4,  We  compute  defined  by 


(9)  Dii=D 


1  + 


l^jxrn.da 


i  =  1,2,3. 


V  r  J 

5.  The  nutrient  transport  equation  will  be 
dc  1  „  _  d^c 


(10)  ^+_L.Ve.v  =  2:D-. 

Ca  7tv  i=i 


5X; 


where:  c  =  c(t,x)  is  the  concentration  of  the  nutrient,  v  =  v(x 
field  of  the  fluid  flow  in  porous  region  (find  by  Chen  et  al.). 


TT 


is  the  porosity. 


)  is  the  velocity 


6.  If  the  solid  inclusions  are  spheres  then  D*  j  =  D22  =  D33 


If  the  convection  is  strong  and  dominate  diffusion  the  above  algorithm  must 
be  modified  as  follows: 

3'.  In  the  unite  cube  Y  we  have  to  solve  the  following  3-D  boundary  value 
problem: 
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-e.  inY^ 

(11)  V^v^=0  inYp 

=  0  onF 

and  are  Y  periodic 

Cj  is  the  unit  vector  in  Yj  direction  j  =  1,2,3. 

4'.  For  X  in  the  hot  zone  we  have  to  compute  the  velocity  field  v(x,y) 
defined  by 

(12)  v(x,y)=-f±.v 

H  5x. 

where  p  =  p(x)  is  the  pressure  field  obtained  by  Chen  et  al.  for  this  region. 

5'.  For  X  in  the  hot  zone  we  have  to  take  the  average  <v>=<v>(x) 
defined  by 

(13)  <v>(x)=-^  Jv(x,y)dy. 

1^1 

Probably  this  will  be  something  very  close  to  the  velocity  field  v(x) 
obtained  by  Chen  et  al,  for  the  region. 

6 .  For  X  in  the  hot  zone  we  have  to  solve  the  following  boundary  value  cell 
problem 

A,Xi-<v>-V^Xi  =  Vi-<Vj  >  inYp 

(14)  h-VyXi=-ni  onF 


Xi  is  Y  periodic 


JXidy  =  0 

where:  <  v  >=<  v  >  (x)  was  obtained  at  5'  “ 

Vj  =  Vi(x,y)  are  the  components  of  v(x,  y)  obtained  at  5' 
Xi  =  Xi  (x,  y)  is  unknown. 

7'.  For  X  in  the  hot  zone  we  have  to  compute  Dy  defined  by 


8'.  The  transport  equation  in  hot  zone  will  be 


(16) 


1  _  35 

H - Vc-  <  V  >=  2] - 


a 


trax,. 


ID,  * 


dx 


jy 
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9'.  If  in  hot  zone  the  pressure  gradient  Vp  is  constant  then  Dy  obtained  at  T 
are  constant  and  in  the  transport  equation  obtained  at  8'  we  can  use  only  the 

symmetric  part  Dy  =~(Dij  +Dji)  of  the  diffusion  tensor.  This  symmetric 

part  is  given  by 

(17)  D-=Dj 


J(x,Vj  +X|V,)dy 


and  the  transport  equation  becomes 


(18)  ^  +  -Lvc.<v>=tiD:,- 


di  7ly 


i=l  J=l 


Bk  d\ 


2.3.  The  flow  in  the  mushy  zone 


Therese  Levy  in  the  paper  [23]  studies  slow  flow  of  an  incompressible 
viscous  fluid  in  an  an  ay  of  a  great  number  of  small  fixed  particles  assuming 
that  he  particle  size  s  and  the  distance  q  between  two  neighbouring 
solids  are  such  that  s  «  q  «  1.  This  assumption  corresponds  to  low 
concentration  of  solid  inclusions.  She  proves  that  there  exists  a  critical  size 
for  particles,  for  which  Brinkman’s  law  occurs.  For  larger  particle  sizes  the 
fluid  flow  is  governed  by  Darcy’s  law  and  smaller  sizes  solids  do  not 
influence  the  flow. 

According  to  Levy’s  paper  if  in  the  mushy  zone  the  size  of  solid  inclusions 
is  subcritical  then  the  steady  flow  in  mushy  zone  is  the  global  convective 
flow  in  cold  zone  defined  by  the  equations: 

(1 9)  (vV)v  =  -Vp“  +  Pr  V"  V  -  Ra  •  Pr-  0  •  e, 

(20)  v-Ve  =  V'0 

(21)  Vv  =  0 

where  the  Prandtl  number  Pr  and  the  Rayleigh  number  Ra  are  defined  by 
Pr  =  Ra  =  p  g(Th  -  Tc)  LVaL-li. 

If  the  size  of  the  solid  inclusions  is  critical  then  the  flow  in  mushy  zone  is  a 
Brinkman  flow  governed  by  a  Brinkman’s  law: 
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(22)  |i  •  A  <  V  >=  Vp°  +  p  •  (p  -H-  <  V  > 

where  (p  is  the  volumic  concentration  of  solid  particles  and  H  is  the 
translation  tensor. 

If  the  size  of  solid  inclusions  in  mushy  zone  is  supper  critical  then  the  flow 
in  mushy  zone  is  governed  by  Darcy’s  law 

(23)  |i-(p-H-<  v>=-Vp°. 

We  note  that  in  subcritical  case  is  not  necessary  to  put  coupling  conditions 
but  in  supercritical  case  is  absolutely  necessary  to  put  this  kind  of  condition 
in  order  to  realize  the  continuity  of  the  velocity  field.  In  this  order  we  can 
use  Beavers-Joseph  postulate  and  Rudraiah  slip  conditions  fsee  1241  1251 

[26]).  ^  I  hi  i. 

Once  global  velocity  field  v  and  pressure  field  p°  are  deteraiined  from 
equations  (19),  (20),  (21),  in  super  critical  case  we  use  the  pressure  field  p® 
to  find,  with  formula  (23),  the  Darcy’s  flow  <  v  > .  With  the  above  Darcy’s 
flow  and  the  global  velocity  field  v  using  Rudraiah  slip  conditions  at  the 
mushy  zone/pure  fluid  interface  we  find  the  small  distance  d  away  from  this 
interface  in  porous  side  where  we  make  the  coupling  between  the  global 
flow  V  and  Darcy’s  flow  <  v  >.  * 

As  following  assume  that  the  mean  flow  in  mushy  zone  is  directed  along  the 
Xi  axis.  This  means  that  <vi>  9^  0  and  <V2>  =  <V3>  =  0. 


2.4  The  transport  equation  of  the  nutrient  in  mushy  zone  when 
convection  and  diffusion  balance  each  other 


Using  homogenization  method  we  can  obtain  the  following  transport 
equation  for  the  solute  in  mushy  zone 

(24) 

^  Tty  i=l  H  5x.5x. 


In  this  equation  c  —  c(t,x)  is  the  solute  concentration,  <Vi>  =  <Vi>(x), 
i  —  1 ,2,3  are  the  components  of  the  fluid  velocity  in  mushy  zone  given  by 
(19)-(21)  or  (22)  or  (23)  and  Dy  are  the  components  of  the  dispersion  tensor. 
It  is  understood  that  equation  (24)  concerning  the  case  when  convection  and 
diffusion  balance  each  other  and  Tty  is  the  porosity  related  to  the  volumic 
concentration  cp  of  solid  particles  by  the  relation  Tty  =  1  -  (p. 
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For  symmetric  solid  inclusions  we  have 

f 

(25)  D,=D  + 

V  |ys|r 

where  the  functions  Xi  are  solutions  of  the  Neumann  problem 

(26)  Ax,  =0inR'\Ys 

(27)  Vx,  •h  =  -nj  onF 

i  =  l,2,3  (see  [53]). 

Steady  state  solutions  of  equation  (24)  are  given  by 
1  ^  dc  d'c 

71,^  ax,  '^dx.dx. 

This  equation  generalizes  the  steady  state  mass  balance  equation 

(29)  ^(vVc)  =  Vc 
Pr 

used  for  the  determination  of  the  solute  field. 

In  equation  (29)  Sc  represent  the  Schmidt  number  Sc  =  p/D  and  Pr  the 
Prandtl  number  Pr  =  p/aL. 

Boundary  condition  for  solving  equation  (28)  can  be  formulate  starting  from 
the  location  of  the  mushy  zone.  We  will  assume  that  the  crystal  is  located  in 
the  region  Qi  =  [-Li/2,  Li/2]x[-L2/2,  L2/2]x[-L3/2,  L3/2],  the  +X3  direction  is 
the  growth  direction  and  the  mushy  zones  are  located  in 
a;  =[-L,  /2,L,  /2]x[-L,  /2,L3  /2]x[L3  /2,(L3  /2)+5]  and 
Q[  =  [-L,  / 2,  L,  / 2]  X  [-L,  /  2,  L3  / 2]  X  [(-  L3  / 2) -  6 -L3  / 2] . 

For  symmetrical  reasons  we  will  consider  only  a  part  of  the  mushy  zone 
and  we  will  assume  that  we  have 

(30) 

'5X3  Pr  ^  ^ 

for  X3  =  L3/2  and  (x,,  X2)€[-Li/2,  Li/2]x[-L2/2,  L2/2] 

which  express  conservation  of  mass  at  the  fluid/crystal  interface. 

On  the  rest  of  the  boundary  of  the  mushy  zone  we  put  a  coupling 

condition  between  the  concentration  field  c  defined  by  (28)  and  the 
concentration  field  c  defined  by  (29)  assuring  the  continuity. 

Once  the  velocity  field  in  mushy  zone  is  determined  we  can  find  the 
concentration  field  c  satisfying  equation  (28)  and  the  above  boundary 
conditions.  This  nutrient  transport  equation  in  the  neighborhood  of  the 
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crystal  is  a  new  more  refined  approach  of  the  real  nutrient  transport  in  the 
region. 


2.5  The  transport  equation  of  the  nutrient  in  mushy  zone  when 
convection  dominates  diffusion 


In  the  case  when  convection  dominates  diffusion  and  the  pressure  field  has 
constant  gradient  in  mushy  zone  we  have  the  following  transport  equation 
for  the  solute: 


5c  1  ^  do 

di  '  dx. 


=  SID„. 


dx.dx. 


(see  [67]).  In  this  case  the  components  of  the  dispersion  tensor  depend  also 
on  the  velocity  field. 

If  the  mean  flow  is  directed  along  the  Xi  axis,  what  was  assumed,  then  Dy  is 
diagonal  with  two  independent  components,  that  are  the  longitudinal  Dl  and 
transversal  dispersivities  Dx;  Dl  =  D,i  and  Dx  =  D22  =  D33.  In  the  case  of 
small  concentration  of  spheroid  solid  inclusions  Mauri  in  [75]  obtained  the 
following  numerical  results: 


Fig.4a.  Longitudinal  diffusivity  in  fixed  beds  of  prolate  spheroids  as  a  function  of  the  eccentricity 
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0-0!? 


Fig.4b.  Transversal  diffijsivity  in  fixed  beds  of  prolate  spheroids  as  a  function  of  the  eccentricity 


Eccentricity 


Fig.Sa.  Longitudinal  diffusivity  in  fixed  beds  of  oblate  spheroids  as  a  function  of  the  eccentricity 
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Fig.5b.  Transversal  difUisivity  in  fixed  beds  of  oblate  spheroids  as  a  fimetion  of  the 


eccentricity 


Also  Mauri  in  [60]  finds  analytically  for  a  dilute  lattice  of  uniform  spheres 
*at  Dt  is  proportional  to  (Pe)^  for  small  Pe  and  is  eight  time  smaller  than 

Therefore  in  our  case  equation  (31)  becomes 

(32) 

^  ^1  dK.  I  c!x^  I 


Steady  state  solutions  of  equation  (32)  are  given  by 

(33)  -i_<  V,  >  =  D,  .  —  +d1 

dx.  ^ 


^x^  [ax"  ax" 


Equation  (33)  generalizes  the  mass  balance  equation  (29)  used  for  the 
determination  of  the  solute  field. 

equation  (28^  boundary  conditions  like  for  the 

Once  velocity  field  in  mushy  zone  is  determined  we  can  find  the 
ZCcLmoI'  --sponding 
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2.6  Stability  and  segregation 


For  obtaining  a  single  crystal  in  a  steady  state  regime,  the  solute  quartz 
concentration  in  the  precrystallization  zone  must  be  independent  on  t,  Xi,  X2. 
An  important  dependence  on  these  variables  may  lead  to  a  dendritic  growth. 
This  means  that  in  the  neighborhood  of  the  seed  the  concentration  c  of  the 
quartz  may  depend  only  on  X3.  Using  this  type  of  function  in  equation  (33) 
we  find  that  c  satisfies 

(34)  Dt-^  =  0. 

ux  2 

The  general  solution  of  equation  (34)  is 

(35)  Co(x3)=K,+K2X5 
where  Ki,  K2  are  constants. 

The  constant  K2  is  determined  by  the  boundary  condition  (30)  and  satisfies 

(36)  K3=^^(l-k)(K,+K2-L3/2) 

Therefore  we  have 


(37)  K3  = 


^(1-k) 

Pr _ 

L3  Pe-Sc/ 


(i-k) 


and 


,(>‘3)=K, 


2  Pr 


We  will  analyse  the  stability  of  a  solute  quartz  distribution  in  mushy  zone 
defined  by  (34).  This  means  that  for  (38)  we  search  perturbations  which  tend 
to  0  for  t  tending  to  +00.  The  reason  is  to  find  time  dependent  regimes  which 
are  able  to  cormect  a  given  dopant  distribution  in  mushy  zone  with  the 
distribution  (38). 

For  the  beginning  we  consider  perturbations  of  the  form 
(39)  Ci(xi.t)  =  Ci(xi>‘^. 


Substituting  (39)  in  (32)  we  find 
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(40)  D 


d^Ci 

^  dx^ 


1 

- <  V,  > - -  -  (TCj  =  0 

Uy  dXj 


and  we  find  that  Co(x3)  given  by  (38)  is  not  stable  with  respect  to  any 
perturbation  of  the  form  (39). 

If  we  impose  to  Ci(xi)  the  condition 

(41)  Ci(-Li/2)  =  c,(L,/2)  =  0 

then  we  obtain  the  following  inequality  for  a: 

(42)  G<- 


<  Vj  >■ 


4(7Cy)^Dl 
From  (40)  and  (41)  we  obtain  also  that 

(43)  c]'(x,,t)=Kr -sinj 


2n7i 

vV 


,a|X, 


•e'""*  n  =  1,2,3,.... 


where  an  is  given  by 
1 


(44)  a„=- 


4D. 


<Vi>-  16DL-7c^-n' 

~ -  + - k - 


L  ("v)= 


'I 


n  =  1,2,3,... 


and  tti  is  given  by 

(45) 


a,  = 


<  Vj  > 


2  •  Tty  •  Dl 

The  concentration  profile  defined  by  (38)  is  asymptotically  stable  with 
respect  to  the  perturbations  (43).  More  than  that,  these  solute  quartz 
distributions  are  asymptotically  stable  with  respect  to  the  perturbations  of 
the  form  Ci  =  Ci(xi,t),  which  satisfy  the  perturbation  equation 


(46) 


1 


dc 

— + 
dt  Tty 


<Vi  > 


5Ci 

5xi 


=  D. 


d^Ci 

dx.f 


and  has  the  property  that  at  the  initial  moment  t  —  0  could  be  expanded  as: 

00 

(47)  c,(x,,0)  =  ^K„ -sin 


Next  we  study  the  stability  of  the  quartz  distribution  given  by  (38)  with 
respect  to  perturbations  of  the  form 

(48)  C2(x2,t)=C2(x2)-e'^. 

In  this  c^e  is  also  stands  that  the  distribution  cq  =  Co(x3)  given  by  (38)  is  not 
stable  with  respect  to  any  perturbation  of  the  form  (48). 
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But  if  we  impose  again  the  spatial  periodicity 
(49)  C2(-L2/2)  =  C2(L2/2)  =  0 

it  turns  that  the  concentration  profiles  Cq  =  Co(x3)  given  by  (38)  are 
asymptotically  stable  with  respect  to  the  perturbations: 

(so)  C2(x2,t)  =  K5 -sin 


2n7t 


where  (?„  are  given  by 


(51) 


4 -It' 


n 


=- 


D, 


The  stability  with  respect  to  the  perturbation  of  the  form 
(52)  C3(x3,t)  =  C3(x3)-e‘^‘ 

is  treated  on  a  similar  manner,  reaching  the  conclusion  that  there  are 
asymptotically  stable  with  respect  to  the  perturbations: 

2n7C 


(53)  c5(x3,t)  =  K3 -sinj 
where  an  are  given  by 

(54) 


(x3-h  +  L3/2) 


and  h  is  the  thickness  of  the  mushy  zone. 

If  we  investigate  the  stability  of  the  solute  quartz  distribution  given  by  (38) 
with  respect  to  the  perturbation  of  the  form 

(55)  c(xi,X2,t)=c,2(xi,X2)-e‘^,  which  satisfies 

(56)  Cj2(“ Lj  /2,X2^—  ^12(^1  / 2, X2)=  Cj2(xj L2  /2)=  Cj2(xj ,L2  /2)=  0 


we  find  that  the  distribution  Cq  =  Co(x3)  given  by  (38)  is  asymptotically  stable 
with  respect  to  the  perturbations: 


(57)  c™(x„Xj,t)  =  Kl"-K; 


'^2m7t  ^ 

^2n7C  ^ 

sin 

L 

J 

•  sin 

T  ^ 

J 

where  a^n  are  given  by 

1 


(58) 


'mn 


4Dr 


<  V|  > 


167t^-m^DJ 


lL  Tty 


47c^n^ 


Dn 


Finally,  if  we  investigate  the  asymptotic  stability  of  the  quartz  distribution 
Co  =  Co(x3)  given  by  (38)  with  respect  to  a  perturbation  of  the  form 

(59)  c(xi,X2,X3,t)=c(xi,X2,X3)-e'^ 
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we  conclude  that  they  are  asymptotically  stable  with  respect  to  the 
perturbation 


(60)  c(xi,X2,X3,t)=K| 


m  ^  n 


Sin 


2in  ] 

'^2m7u 

L 

V  ^1  J 

•sin 

L 

^  ^2  J 

sin 

where  ciimn  are  defined  by: 

1 


2n7c 


-(x3  -  h  +  L3  /2) 


(61) 


^Imn 


4D, 


<Vj  >' 
TCv 


1671^^^dJ 


47i^m^ 


2 

^2 


■  Dt  - 


47i"n' 

L\ 


•Dn 


The  stability  results  obtained  till  now,  could  create  the  impression  that  the 
distribution  Co(x3)  given  by  (38)  is  stable  with  respect  to  any  perturbation. 
For  the  reason  for  which  we  make  this  study  this  would  be  good  but 
unfortunately  this  is  not  true.  We  showed  that,  generally  speaking,  for  any 
a  >  0  and  8  >  0  there  exists  a  perturbation  Ci(xi,t)  =  Ci(xi)e*^  such  that  for 
t  —  0  we  have  supjc,  (x,0)|  <  £■ .  In  fact,  this  was  the  reason  why  we  have 

reduced  the  family  of  the  perturbation  of  the  form  Ci(xi,t)  =  Ci(xi)e^  and  we 

have  considered  only  perturbations  for  which  we  have 

Ci(-Li/2)  =  c,(Li/2)  =  0.  It  must  be  noted  that  the  vanishing  condition  of 

Cl  =  Ci(xi)  in  -Li/2  and  Li/2  could  not  be  substituted  only  by  a  periodicity 

condition. 

In  the  problem  of  the  stability  of  the  solute  quartz  distribution  Cq  =  Co(x3) 
with  respect  to  perturbations  of  the  form  (55),  the  vanishing  condition  on  the 
margins  |x2|  =  L2/2  (ci2(xi,-L2/2)=Cj2(xj,L2/2)=0)  can  be  replaced 
with  a  periodicity  condition  Ci2(x,,-L2/2)  =  Ci2(x,,L2/2),  but  we  must 
conserve  the  vanishing  condition  on  the  margins  ' 
1^1 1  =  L2/2  (Ci2(-  L2  /2,X2)=  Ci2(L2  /2,X2)=  0). 

In  a  similar  way  we  can  prove  that  the  quartz  distribution  cq  =  Co(x3)  given 
by  equation  (38)  is  asymptotically  stable  with  respect  to  the  following 
perturbations: 


(62)  ^n,A,,Y (^1 5 X2 j t)— Kj  •  sin 


2n7t 


sh 


Li 

^2(^2 


1 

V  y 

-1 


le“'^>  .y 


wheren=  1,2,3, ...,yeR',y;6  0,A,€R’;;.<0;  Kl*  gR'; 
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A,>-- 


4D, 


<Vi>^  16%^dI 

+  - - 


'I 


;zi.2=+ «i  = 


4D, 


<v,  1671^ -n^-D? 


+  • 


Tt^ 


D, 


-X 


<Vj  > 

iTty  -Dl 


(63)  Cn^(x,,X2,t)=K{’ -sin 


2mt 


Xj  e“'^' 


where  n=  1,2,3,...,  yeR',  y?i  0;  K"  eR';  a,  = 


<  Vj  > 


= 


4D, 


<v,  1671^  -n^  -  Dl 


+  ■ 


7i; 


(64)  c„;,.y(xi,X2,t)=Kf -sin 


2n7t 


V  ^-1  J 


e^i'’^'  .y. 


P2L2 


COS 

2 


CTn  51  t 


■COSP2X2  'S  "’*■ 


where  n  =  1,2,3,...,  yeR',  y^O,  A,eR';  A.  >  0;  K"  eR'; 


X^^Dj;  P2=. 


A 

'Dt 


a,  = 


<  Vj  > 

27tYDL 


^n,A. 


4D, 


<Vi  ^  167C^  -n^  -D^ 


A 


-X. 


The  perturbation  (62)  corresponding  to  the  value  Xc  given  by 

2 


(65)  >.,=■ 


1 


4D, 


<  Vj  > 


TC. 


167t^ 


is  critical  in  the  sense  that  for  the  values  of  X  smaller  than  Xc  the  distribution 
(38)  of  the  solute  quartz  is  not  stable  with  respect  to  the  perturbations  (62). 
Also  in  the  case  of  perturbations  of  the  form  c(xi,  X2,  X3,  t)  =  c(xi,  X2,  X3)-e®* 
the  vanishing  conditions  on  the  boundaries  [xj]  =  1^/2  and  X3  =  h  -  L3/2,  X3=h 

can  be  replaced  by  periodicity  conditions,  maintaining  the  vanishing 
conditions  on  the  boundary  |x,|  =  Z,  /  2 . 


47 


2.7  Conclusions 

'■  identified  transport  equation,  which  generalizes  the  convective- 

di^sive  mass  transport  equation  used  in  general  in  crystal  growth  for  the 
determination  of  the  solute  field. 

li"  ®‘l“^*ions  describe  the  dispersion  of  the 

solute  when  convection  and  diffusion  balance  each  other  respective  when 
convection  dominates  diffusion.  These  equations  take  into"ac“e 
Structure  which  already  exists  in  mushy  zone,  the  interaction  between 

complicated  global  How  patterns  and  the  convective-diffusive  dispersion 
mechanism  in  this  region.  ^i^persion 

"  M°ASi”p.::mprz^ 

5.  It  must  be  noted  that  these  equations  do  not  take  into  account  the 

m  erac  ion  of  the  solute  field  with  the  crystal.  From  this  point  of  view 
these  equations  can  be  improved. 

6.  The  stability  calculus  suggests  that  there  are  a  lot  of  time  dependent 
regimes,  which  tend  to  a  steady  state  uniform  distributed  concentration 

^  boundary  control  calculus,  presented  for  Bridaman 

Stoebarger  growth,  suggest  in  how  way  we  can  reduced  the  segregation 
m  hydrothermal  growth  system  also.  segregation 
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